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Need for automated methods 


for image registration 
Large data 
sets 


Launch of several 
planetary missions 
Design of new and 
powerful sensors 


Objective 


¢ Crater detection in planetary images 
¢ Development of an image registration method based on the 


extracted features 
Crater Crater 
Detection Map 
Crater 
Map 
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Crater Detection 


Marked Point Process Model 
Energy Function 

— Multiple Birth and Death Algorithm 
— Region-of-Interest Approach 
Experimental Results 


Marked Point Processes 
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Crater detection based on a marked point process (MPP) model 


Realizations Configurations of objects, each 


MPP: Stochastic Process described by a marked point 


Mathematical Formulation 


A point process X, defined over a bounded subset P of R* maps from a 
probability space to a configuration of points in P. 


Realizations of the process X are random configurations x of points, x = {xj,..., xy}, 
where x; is the location of the i‘ point in the image plane (x; € P) 


A configuration of an MPP consists of a point process whose points are enriched with 
additional parameters, called marks and aimed at parameterizing objects linked to the points. 


Bayesian approach: Maximum a posteriori (MAP) rule to fit the model to the image is equivalent 
to minimizing an energy function (computationally challenging) 


Marked Point Process for Crater Detection fy 


m Modelling 
<—= 


Distribution of craters Distribution of ellipses 


on the surface on the image plane 


4 Realization Example 
Center: (0,¥0) Point 
Point Process 
Major Axis: a RES TX iy wig ha} 
Craters 
Minor Axis: b Marks 


Parameters 


Orientation: @ Xi = (Xi Voi A, Dj, 9;) 


Crater Detection — Energy Function 
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Energy function of the configuration X = {x;, x2,...,X,} wrt the 
extracted set C of contour pixels (Canny): 


U(X|C) = Up(X) + U,(C|X) 


Prior Likelihood 
Repulsion coefficient based on the Two terms, one based on a correlation measure, the other based on 
overlapping of the ellipses (overlapping a distance measure (fit between contours and realization of X) 
craters are quite unlikely) 
n 
dyp(x?,C) |xPnc| 

1 Xx; AX; U,(C\X) = > CHM — |e 
U(Xy=- ) a “arr (cl 

n x; V xj i=1 

xjAxj>0 


x? = set of pixels corresponding to ellipse x; in the image plane 
x; V x;= area of union of ellipses x; and x; | d3¢(x?,C) = Hausdorff distance between ellipse x; and the contours: 
x; A x;= area of intersection of x; and x; 


dy,(A, B) = ees hel inf d(a, B); oe je d(a, BD 


acA BEB 


Classical distance between sets d(A, B) = 0 


Crater Detection — Energy Minimization 
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Markov chain sampled by a 
multiple birth and death 
(MBD) algorithm 


Markov chain Monte Carlo-type method > 
Simulated Annealing scheme 


a Contour Extraction (Canny) and Parameter Initialization 
Initialization 


Generation of the Birth Map used to speed up the convergence 


S Computation of the Birth Probabilities for each pixel 
Birth Step 
Birth Phase 


Energy Computation for all the ellipses 
Death Step Computation of Death Probability and Death Phase according to Likelihood 
Configuration refinement according to Prior 


Convergence Test 


Update + 


MBD - Birth and Death Steps 


Goddard Space 
Flight Center 


Birth Step Death Step 
For each pixel s in the image, compute the For each ellipse x; in the configuration, 
birth probability as min{6é - B(s), 1}, where: compute the death probability as d(x;): 
b(s) 6 - a(x;) 
Se d(x;) = —————_-~ 
B(s) ¥, b(s) (xi) 1+6-a(x;) 


b(s) is the birth map computed from the a(x;) = exp|—B(U, (X\{x}1C) — U_(X|C))| 
contour map using generalized Hough transform 
and Gaussian filtering 


NASA Crater Detection — Region Based Approach i 


¢ MBD is computationally heavy 


Region-Based Approach ¢ Computational burden increases 
with image size 


Region Based Flowchart and Example 


Initialization 


Crater Detection — Data Sets 


Goddard Space 
Flight Center 


¢ 6 THEMIS (Thermal Emission Imaging System) 
images, TIR, 100m resolution, Mars Odissey mission 


¢ 7 HRSC (High Resolution Stereo Color) images, VIS, 
~20m resolution, Mars Express mission 


¢ Image sizes from 1581 x 1827 to 2950 x 5742 pixels 


Quantitative Performance Assessment of the crater detection algorithm: 
Detection Percentage (D), Branching Factor (B), and Quality Percentage (Q) 


TP FP TP 
Data p= i; ee 0 = —______ 
TP +FN TP TP+FP+FN 
Avg on all THEMIS 0.91 0.10 0.83 
Avg on all HRSC 0.89 0.06 0.85 
Avg on all images 0.90 0.09 0.84 


> 


HRSC Sensor 


HRSC Sensor 
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Crater geometric properties extracted by the proposed method 


Crater 
Crater 1 
Crater 2 
Crater 3 
Crater 4 


Crater 5 


C = (Xo, Yo) 
(139, 393) 
(258, 756) 
(343, 23) 
(591, 215) 
(919, 157) 


Semi-axes (a, b) 
(35, 33) 
(51,50) 
(13,12) 

(19, 18) 
(15,14) 


THEMIS Sensor 


Orientation 0 
64° 
Ts? 
180° 
31° 
106° 


Image Registration 


— 2-step Approach 
— Experimental Results 


Crater-based 
Registration 


Feature-based registration 

¢ Min Hausdorff distance (d3,) between extracted 
craters through genetic algorithm 

¢ Fast but sensitive to accuracy of crater maps 


Area-based registration 

¢ Max Mutual Information (MI) through genetic 
algorithm 

¢ Highly accurate but computationally heavy 


Image Registration — 2-Step Optimization 
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RST (Rotation-Scale-Translation) 
transforms p = (ty, ty, 0, k) 


Ml-based 
Registration 


Refinement: registration based on 
mutual information ina 
neighborhood of p > p* 
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Transformation found for an interactively 
selected region of interest > pz 
Dp u (ty, ty, 0, k) ae —k cos(@) x» —k sin(@) yo + ty + Xp 
cca ae k sin(@) x9 — k cos(9) yo + ty + Vo 
0 


Transformation derived for the entire k 


Image > pa 
Pa = Cie Ty, B, a) 


DEO WES = Superposition of Reference and Input 


Input image 


Image Registration — Data Sets fy 


Semi-simulated image pairs Real multi-temporal 


. image pairs 
20 pairs composed of one real 


THEMIS or HRSC image and of an Real multi-temporal 

image obtained by applying a pair of LROC (Lunar 

synthetic transform and AWGN Reconnaissance 
Orbiter Camera) 


Quantitative validation with respect to images 


the true transform (RMSE) 
100m resolution 


Source Image 


Only qualitative visual 
analysis is available, 
as no ground truth is 
available 


Registration 
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RGB of original non- 


registered data 


1396x2334 HRSC 


Data set RMSE [pixel] Der 
THEMIS (10 data sets) 0.31 p* 
HRSC (10 data sets) 0.22 RMSE 15 Step 
Average (20 data sets) 0.26 RMSE 2"9 Step 


18 


~~ ___— 


| 


RGB of registered data 


Left Image Right Image 

(7.05, 35.91, 0.18°, 1.071) (76.59, 19.96, 2.17°, 1.031) 

(7.04, 35.92,0.19°, 1.071) (76.41, 20.06, 2.18°, 1.031) 
0.79 0.51 


0.16 0.33 


Registration Results with Real Data 


Visually accurate matching between 
reference and registered images in 
the real multitemporal data set 


Checkerboard 
representation of 
the registered 
images (zoom on 
details) 


Registration Results with Real Data 


Visually accurate matching between 
reference and registered images in 
the real multitemporal data set 
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Conclusion 


Conclusions 


Accurate crater maps, useful for both 
image registration and planetary science, 
were obtained from data from different 
sensors. 


Higher accuracy as compared to 
previous work on crater detection (not 
shown for brevity) 


Reduced time for convergence thanks to 
a region-based approach 


Sub-pixel accuracy and __ visual 
precision in registration: effectiveness of 
the proposed 2-step registration method 


Conclusions and Future Developments 
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Future Developments 


Test in conjunction with a_ parallel 
implementation (e.g. computer cluster) 


Validation with multi-sensor real 
images 


Extension to other applications 
requiring the extraction of ellipsoidal or 
circular features, e.g. optical Earth 
observation images or medical images 
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For each pixel in the image compute the Birth Probability as min{6 - B(s), 1}, where: 


b(s) 
eae) 


Being b(s) the Birth Map computed from the Canny Contour Map 


MBD — Death Step 
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For each ellipse x; in the configuration compute the Death Probability as d(x;), where 


6: a(x;) 


aC 1+6:-a(x%;) 


and a(x;) = 2 B(UL(AxiMg)-ULUg) = o8-Ui(zillg) 


The complete Flowchart of the Death Step is as follows: 


Hausdorff Distance 


N& P 


Similarity = mean, > > [anes xr) 


i=1 t=1 


c = craters in Input Image 

N° = sum(pixels in crater c in Input Image) 

P = sum(craters’ border pixels in Ref Image) 
x; = coord of pixel iin crater c in Input Image 
X, = coord of pixel t in Ref Image’s craters 


MI(X,Y) = ». 


Similarity Measures 


Mutual Information 


Pxy (x,y) 
Dx (x) py(y) 


». Dxy(x,y) log ( 


xEX yEY 


X: pixel intensity in Reference Image 

Y: pixel intensity in Input Image 

Px (x): probability density function (pdf) of X 
Py (y): probability density function (pdf) of Y 
Px y (x,y): joint pdf of X and Y 


Dx (x) Py) Pxy (x,y) 


estimated through 
the corresponding 
image histograms 
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Rotation — Scale — Translation Transformation 


Transformation vector 


p= (ty, ty, 0, k) Matrix Formulation 
ft,,ty}:T lati ; d kcos(@) ksin(@) t, ‘ 
oily. ranslations in x an y Gy) = Enc). Iecasco) @ y 
0: Rotation angle y/ 11 


k: Scaling Factor 


Original Rotation Scaling Translation 


Region of Interest Approach fy 
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xX I,(X,Y), Ip(x, y): Two Images 
Iz: sub-image of I, such that Ip(0,0) = I4(%o, Vo) 


Da = (Ty, Ty, B, a): RST transformation vector transforming I, into re 
Dp = (ty, ty,0,k): RST transformation vector transforming I, into [5” 


I (0,0) = Ve (Xo, Vo) 


transformed 


; ; Transformation: pp a te eee 
y Given: m ee oe ee (Cen) Find: Transformation: pa 
F 


Expressing the transformation in Matrix Form 


From the image acos(B) asin(B) T, ro . 
Cay T,, = & ee Ty (X,Y) = (X4Y%) This should also hold 
¥ Tp ,(* + Xo, ¥ + Yo) = (X' + Xo, ¥' + Yo) 
P=yT Vo _(kcos(@) ksin(@) ty\ eee - 
"be =\_k sin(@) kcos(@) ty) ° Typ) = Y) 
Plugging T,,, into this equation and replacing a 
x' and y’ according to Ty, Knowing a = k and solving in P, = (0,0) and P, = (—x9, —yo) 
k cos(@)x +k sin(@) y+ ty +x = 
é —k cos(@) xp — k sin(@) yo + ty + Xp 
a cos(B)(x + x9) + a sin(B)(y + yo) + Ty . 
B ° BY + Yo ne k sin(@) x9 — k cos(@) yo + ty + Yo 
—k sin(9) x + k cos(@) y+ ty + Yo = : 


a sin(B)(x + x9) + a cos(B)(y + Yo) + Ty 


RMS Error Computation fy 
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Erorr Transformation 


Pe = (tre, byexOoy ka) ? Qp, = Qn : Gees 


k 
—» ke=7 0 = 2-H 
1 


tre = ty — ke(ty1 cos(O,) + ty: sin(62)) 
tye = ty2 — ke(ty1 cos(@,) — ty: sin(9.)) 


Ground Truth Transformation 
Pet = (tas ty1,94, k,) BD i) = Qa * ii 


ComputedTransformation 
p= (ty, bia @, k) 2 130%, 9) = 05" (%¥, 1]7 


(x,y) € Image, [x’, y’, 1]” = Qp,: [x,y, 1]" — | ite eee er 5 " [ee 


RMS Error: E(p,) = \— J mea —x)2 + (y’ — y)2dxdy, a = A* + B? 


1 A -B 
E*(p.) = al | (k, cos(@,)x + ke sin(0,) y + tye — x)* + (—k, sin(O,) x + k, cos(@,) y + tye — y) dxdy 
o Jo 


a 
E*(p.) = 5 (k2 — 2k, cos(@,) + 1) + (t2, + t2,) — (At2, + Bt2,)(1 — ke cos(Oe)) — ke(Atye — Btxe) sin(e) 


